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The spin-boson model is a paradigm for studying decoherence, relaxation, entanglement and 
other effects that arise in a quantum system coupled to environmental degrees of freedom. At 
zero temperature, a localization-delocalization phase transition is known to exist in the sub-Ohmic 
regime, where the standard density matrix renormalization group algorithm is inadequate due to the 
divergence in the number of low-frequency modes. This limitation is circumvented in this work by 
symmetrically optimizing the phonon basis and introducing an order parameter accounting for the 
U(l) symmetry for a two-bath spin-boson model, by which we are able to determine the classification 
and criticality of the phase transition explicitly. Compared with variational results, the critical phase 
is characterized by spontaneous vanishing of boson displacements in both the baths, resulting in an 
accurate phase diagram with three model parameters. 

PACS numbers: 


Introduction- Much attention has been devoted in re¬ 
cent years to optical properties of natural photosynthetic 
systems and organic photovoltaic devices [M3, 

where the quantum aspect of excitons and phonons is 
increasingly recognized to be essential in boosting the 
power conversion efficiency. In organic systems, e.g., de- 
localization of wave functions is found to be essential for 
the dissociation of excitons Q. As an excitonic paradigm, 
a two-level system is described by a single spin one-half 
which is coupled to its phonon environment represented 
by boson modes. This leads to the celebrated spin-boson 
model (SBM) 0 0 , providing avenues to study the 
phase transition between the localized and delocalized 
phases, or from a slightly different perspective, the dy¬ 
namical phase transition between the coherent and inco¬ 
herent phases [Tdl - lzi . Despite its simplicity, the SBM is 
a highly nontrivial model in nearly all aspects, and cur¬ 
rently contention still surrounds the existence and the 
precise locations of phase transitions. Designed to help 
understand intrinsic mechanisms of coherent exciton dy¬ 
namics, every theoretical approach typically works ac¬ 
curately only in a certain applicable regime, which is far 
away from the critical point, preventing the method from 
addressing issues related to the phase transition. 

The density matrix renormalization group (DMRG) is 
a powerful numerical technique to study the low-l yin g 
states in strongly coupled, one-dimensional systems [231] . 
Similar to that for the numerical renormalization group 
f2ft method, the orthogonal polynomials theory can be 
employed to map the SBM to a one-dimensional chain 
with only nearest-neighbor coupling [0 , allowing one to 
straightforwardly adapt the DMRG method to study the 
SBM with only model-free approximations. Over the last 
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few years, this approach has been extensively used for de¬ 
tailed studies of the phase transition of SBM i26l43TI |. Our 
recent work was devoted to the two-bath SBM (TBSBM) 
to investigate the phase transition in a comprehensive 
manner [30l 1311] . We have also examined the dynamics of 
SBM with the time-dependent DMRG (t-DMRG) algo¬ 
rithm, compared it with two other established methods, 
and demonstrated that a unitary transformation for the 
state yields reliable, accurate results [32l. l33j. 

The approach of the optimal phonon basis, originally 
developed to deal with coupled electron-phonon systems 
[34[, is often adopted to reduce the dimension of the Fock 
space of bosons (28|, 30] ■ Although it has been utilized for 
a variety of models l3M3, a serious problem arises in 
the context of the TBSBM wherein the symmetry is nu¬ 
merically broken and the analysis of the phase transition 
becomes obstructed [0. In this work, we circumvent this 
problem by symmetrically optimizing the phonon basis 
and constructing numerics-friendly operators. Our ap¬ 
proach is adapted specifically for the TBSBM such that 
the doubly degenerate ground states can be obtained in 
a credible manner, and properties of the phase transition 
can be studied with sufficient precision. 

Model and methodology- We consider the TBSBM 
in which a single spin is coupled diagonally and off- 
diagonally to two independent baths characterized by 
continuum spectral densities. The corresponding Hamil¬ 
tonian can be written as 


* = EE 

is=z,x l 


uib\ v bi, v + — A i^{b\ v + 


(1) 


where a z and a x are the Pauli operators, b\ u [bi tV ) is the 
creation (annihilation) operator of the Z-th mode of fre¬ 
quency oji in the v-th bath (v = z, x), and A i t „ represents 
the corresponding spin-bath coupling strength. In the 
traditional SBM, the bath spectral density has a cut-off 
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frequency ui c . For simplicity, the same cut-off frequency 
is assigned to the spectral density functions of the two 
baths in this work, i.e., J„(w) = 2na u ujl~ s io s e~ ul ^ Uc with 
a v being the dimensionless coupling strength for the v 
bath and s being the exponent. The case of s < 1 corre¬ 
sponds to the sub-Ohmic regime in which the localized, 
critical and delocalized phases have been claimed to exist 
[28j . We will focus on the sub-Ohmic regime due to its 
relevance and complexity. 

As discussed earlier, the symmetry in the Hamiltonian 
CD is of paramount importance to the numerical preci¬ 
sion. To facilitate discussion, we introduce the operators 


Of = ±cr 2 e” Ei , Of = ±a x e in Ei b L b ‘^, ( 2 ) 


which commute with the Hamiltonian. Also of interest 
is their product, i.e., 


q C _ OtOf. = i( a y e™^ b U bl ’* +b l* b, ’* ) , 


( 3 ) 


with £, tj>, if = ± following the common rule of prod¬ 
ucts. Together with the identity operator it can 
then be verified straightforwardly that the eight opera¬ 
tors X ± , Of , Of , Of form a non-abelian group G, and its 
center is represented by {2A}. The factor group G/fl^} 
is an abelian group whose irreducible representations are 
given by four one-dimensional ones, indicating the U(l) 
symmetry when a z = a x [29j ] . On the other hand, the 
two-dimensional representation of the non-abelian group 
G, characterized by a nontrivial central extension of its 
factor group, participates in the decomposition in irre¬ 
ducible re pre sentations, resulting in the Z 2 symmetry if 
a z f a x |3l|. Subsequently, eigenstates of the system 
and, the ground state in particular, are doubly degen¬ 
erate, a novel feature that allows specifications of the 
numerical precision in dealing with the symmetry. 

We thus proceed to develop the DMRG algorithm 
to deal with the highly symmetrical model. Widely 
used to study the SBM and related models I2t3432l[. the 
DMRG approach starts with the discretization of the bo¬ 
son modes, and employs the orthogonal polynomials the¬ 
ory to represent the renormalized modes by a set of boson 
sites p5l |. with a transformed Hamiltonian 


* = E 


v=z.x 
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AOPB SOPB 


FIG. 1: Schematic for the two kinds of optimal phonon basis, 
(a) The AOPB is shown which takes into account the phonon 
bases with shifted displacements f|s|. (b) The SOPB is shown 
which considers the globally lowest phonon bases. 


Fock space and retain a finite number of bare phonon 
states for each mode, resulting in a so-called restricted 
phonon basis. Previous tests [13 reveal that the re¬ 
stricted basis method is applicable for the SBM away 
from the critical point, but fails to capture the phase- 
transition properties in the vicinity of the critical point. 
An approach employing an optimal phonon basis (OPB) 
was adopted by Zhang et al., yielding improved results 
341. Details of the OPB-adapted DMRG algorithm can 
be found in the Supplementary Material. Numerical dif¬ 
ficulties arise, however, in the OPB approach while cor¬ 
rectly tackling the symmetry issues. A seemingly simple 
solution is to add an infinitesimal bias and to approach 
the critical point asymptotically. However, this would 
lead to a disastrously large difference in symmetry be¬ 
tween a z a x and a z = a x (the critical point) as dis¬ 
cussed earlier. To circumvent this dilemma, we propose 
two techniques to enforce the model symmetry. 

It is realized that in the implementation of the OPB 
approach, many low probability states are discarded. As 
sketched in Fig. Oja), e.g., if the sign of the calculated bo¬ 
son displacement is negative, those states with a positive 
sign must be eliminated as they are almost orthogonal 
to the negative-displacement states. To recover the sym¬ 
metry, therefore, the information in the eliminated states 
must be fed into the reduced density matrix, as depicted 
in Fig. [Ub). In particular, if |gf is the calculated ground 
state with i as the index of the left free site, we apply 
the parity operator Vi(= e* 7rb ‘ 6< ) onto \gi), obtaining \g[) 
with an opposite sign of the boson displacement on i-tii 
site. This is numerically feasible as the i-th site is the 
one that has the Fock space based upon the bare phonon 
basis. With the two states obtained, the reduced density 
matrix of the i-th site is calculated by 


Pi = Tr E [a\gi)(gi\ + (1 - a)\g\) {g'fW 


( 5 ) 


where rp, is the renormalized coupling calculated from 
r/v = f 0 ° Jv(aj)dui. Herein, 04 and ti are the frequency 
and the hopping integral for the *-th site of bosons, re¬ 
spectively, and the expressions for them can be found in 
Refs. [iHij]- 

Despite discretization of the spectral densities, the 
number of the bare phonon basis in the local Fock space 
for each renormalized boson mode is still infinite, hinder¬ 
ing numerical calculations. A remedy is to truncate the 


where Tr E denotes the partial trace over all the sites ex¬ 
cept the i-th one and a is the portion of the state | gf in 
the reduced density matrix. In the absence of bias, it is 
intuitive to set a to 0.5. We then name the new basis ob¬ 
tained from the adapted reduced density matrix as the 
“symmetrically optimized phonon basis (SOPB),” and 
accordingly, the conventional OPB is called the “asym¬ 
metrically optimized phonon basis (AOPB).” We will 
show that the optimization procedure based on the SOPB 
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yields more accurate results. 

With one run of the DMRG algorithm, one degener¬ 
ate ground state is obtained, and with a second run, 
we may obtain another, allowing for a linear combina¬ 
tion of the two [29|]. However, symmetry breaking still 
persists reducing the accuracy and the efficiency. A so¬ 
phisticated approach would then be to apply the oper¬ 
ators © to the calculated ground state. In order to 
implement this operation, through the parity operator 
V = V t = e' i7r 5A b i h *, a unitary transformation on all 
bosonic modes must be applied to flip the sign of dis¬ 
placements. As the bases of the boson modes have been 
symmetrically optimized, however, the number operator 
hi = b\bi ceases to be diagonal. Moreover, the two de¬ 
generate states are almost orthogonal with each other, so 
that a simple diagonalization for hi can no longer guar¬ 
antee numerical precision. In this context, we introduce 
a more numerically friendly treatment of the parity op¬ 
erator by recasting it in the form 

Vi= ^2 Wi){m\+ e * 59 | n i)( n *l> (6) 

m=even nj=odd 

where | rii) is the eigenstate of the number operator hi at 
i-th site, and 59 is a small angle. Herein, following the 
t-DMRG algorithm [38], the angle n is divided into many 
small steps (59) and the operator is applied incrementally 
to the ground state. The operator does not act on the 
even-number states, while for the odd-number states one 
can make Vi(59) act cumulatively onto the state until a 
certain angle 9. If 9 = n the action is equivalent to that 
of the parity operator, and if 9 = 27t it is an identity 
operator. With this approach, reliable results can be 
obtained for all model parameters. 

Results- In all the calculations we have carried out, 
the number of transformed bosonic site is set to be 50, 
the number of bare phonon basis is 16, and the DMRG 
truncating number is 64. Within these parameters, the 
error is reduced below 10 -5 . The computation is time- 
consuming, e.g., on a single 2.13 GHz processor one run 
for a set of model parameters needs more than one hun¬ 
dred hours of CPU time. 

We first discuss the deep sub-Ohmic regime with s < 
0.5, in which a transition from localized to delocalized 
phases has been discussed in our previous work 30]. 
Fig. 0[a) shows |(ct 2 }| and the ground-state energy E g 
for various values of a x with s = 0.25 and a z = 0.02. We 
compare three cases, with SOPB, with AOPB, and with¬ 
out any OPB adaption, and present results in the vicinity 
of the critical point. With |(cr 2 )| plotted as a function of 
the a x for the three cases and compared to our previous 
work S3, 53 , a much sharper decrease of |(cr 2 }| from a 
finite value to zero is found after the implementation of 
the parity symmetry. An a x increment of 0.0002 is taken 
around the critical point, which is numerically equivalent 
to being infinitesimal. If the approach with SOPB gives 
rise to more precise results, it is implied that the phase 
transition here is more likely to be of first order. In ad¬ 
dition, the ground-state energy is shown in the inset for 



FIG. 2: |(<r 2 }| versus a x for three cases of the basis with 
(a) s = 0.25, a z = 0.02 and (b) s = 0.6, a z = 0.1. The 
green dash-dot lines indicate the phase boundary at a x — a z . 
The insets show the ground-state energy E g versus a x . The 
dotted lines in the inset of (a) are guides for eyes of the linear 
relationship. 


the case with SOPB with the parity symmetry fully con¬ 
sidered. It is observed that there is an obvious kink at 
the critical point implying the phase transition. To check 
the precision of the SOPB-adapted method, we show in 
the Supplementary Material the phase-angle dependence 
of the bosonic displacement, which serves as a measure of 
the precision. The results readily demonstrate the signif¬ 
icant benefits accrued from the precision improvement. 

The shallow sub-Ohmic regime with 0.5 < s < 1 
presents even richer physics. It has been claimed that 
when s > 0.75, there is a so-called critical phase at 
q 2 = a x with (cx 2 ) and (a x ) spontaneously vanishing 
[28j- From a mean-field analysis (3lj| . however, a similar 
phenomenon is found for s > 0.5 instead of s > 0.75. The 
discrepancy may be attributed to the fact that the mean- 
field theory is valid in the weak-coupling limit, while pre¬ 
vious DMRG calculations are applicable in the relatively 
strong coupling regime. To resolve the problem, it is nec¬ 
essary to work with a greater parameter space. 

To this end, we first apply the SOPB-adapted approach 
to the case of s = 0.6. In Fig. [2Kb), we display |(cr 2 )| as 
a function of a x obtained with SOPB, with AOPB and 
without OPB. It is found for the cases with AOPB and 
without OPB, | (cr 2 ) | shows a rather sudden change across 
the critical point. While adopting the SOPB method, the 
curve of the |(cr 2 )| across the recognized critical point 
becomes much smoother as compared to the other two 
cases. This effect implies that the transition is of higher 
order than that with the deep sub-Ohmic bath. More¬ 
over, the curve of the ground-state energy shown in the 
inset ofFig.db) is also smooth close to the critical point. 

In practice, merely considering the magnetization is 
insufficient to determine features of the phase transition. 
Derivation of a more explicit quantity with the complete 
information of both spin and baths is required for this 
purpose. The argument from the group theory shows 
that the operators © are the generators of the parity 
symmetry, with eigen-values of +1 or —1. Following the 
group-theory argument, hereafter we calculate two quan¬ 
tities, (O z ) and ( O x ), involving predictions of both the 
spin and the boson components. The calculation becomes 
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possible because, as tested, the action of the operator V 
is precise based on the SOPB method. 

For further clarity, we define an order parameter as 
C = yj(Oz) 2 + (O x ) 2 - In the localized and delocalized 
phases with a z ^ a x , due to the orthogonality of the 
states O x \g ) and O z \g ), either ( O z ) or ( O x ) should be 1, 
such that the quantity £ will always be unity. In the crit¬ 
ical phase, as stated, the significance phenomenon is the 
spontaneous vanishing of (a z ) and (a x ), which obviously 
results in the vanishing displacements of bosonic modes. 
Let X and Z be the displacements of the bosonic mode in 
the two baths, respectively, which form a X — Z plane for 
all the modes. The critical phase then refers to the case 
in which all the modes in the ground state are located at 
the origin of the X — Z plane, indicating ( O z ), ( O x } and £ 
are all unity. Out of the critical phase, however, the con¬ 
tinuous U(l) symmetry dominates the critical point with 
a z = ol x . It implies that the bosonic modes have certain 
displacements and deviate from the origin, leading to £ 
smaller than one. The order parameter £ accounting for 
the U(l) symmetry is thus a measurement of the effective 
distance of this deviation. 

Fig-El^a) shows the a x dependence of ( with a z = 0.02 
and 0.1 for s = 0.25 and 0.6, respectively. It is clear 
that for the deep sub-Ohmic case, ( undergoes a sud¬ 
den change between 1 and 0 close to the critical point 
(a x = a z ). In terms of the above arguments, this find¬ 
ing proves the absence of a critical phase and clarifies 
the first-order phase transition for the sub-Ohmic bath. 
More interesting is the shallow sub-Ohmic case, in which 
C is between 0 and 1 at the critical point implying the 
possible appearance of the critical phase. In order to see 
the transition point between the localized and critical 
phases, we show in Fig. [3jb) the s-dependent £ for five 
a z 's right at a x = a z . It is found that £ vanishes if s is 
small, while it becomes unity if s is larger than a certain 
value. More importantly, the larger the a z , the larger 
the value of s for which £ equals to one. In particular, 
for a z = 0.1, the transition point ( becomes unity is lo- 


0.05 0.10 0.15 



FIG. 3: (a) The order parameter £ versus a x changing from 
0.015 to 0.035 for s = 0.25, a z = 0.02 and changing 0.05 to 
0.15 for s = 0.6, a z = 0.1. The dash-dot line denotes the 
phase boundary, (b) ( as a function of s for five sets of a z 
and a x . The inset shows the phase diagram consisting of 
the localized phase (LP), the delocalized phase (DP) and the 
critical phase (CP). 



FIG. 4: The order parameter £ calculated by the variational 
approach with rotational optimization shown in the X — Z 
plane for (a) s = 0.25, a z = a x = 0.02 and (b) s = 0.6, a z = 
a x = 0.1. The radius of the circle in the X — Z plane is 
determined by averaged boson displacements. The schematics 
in the X — Z plane shows the displacements with respect to 
the calculated ground states, which are pointed out by the 
arrows as well. 


cated at s = 0.75 which is in agreement with that of the 
previous work [28j . On the other hand, a scaling analy¬ 
sis of a z to the weak-coupling limit shows the transition 
point of £ = 1 moving to s = 0.5 consistent with the 
expectation of the mean-held analysis (3jJ] . Based on the 
transition points of £, we draw a phase diagram shown in 
the inset of Fig. [3j where the boundaries of the localized, 
delocalized and critical phases are explicitly determined. 


Finally, to ensure that the vanishing ( is not caused by 
numerical instability, we have also carried out the calcu¬ 
lations based on variational theory 31| . Via this method, 
we first produce a state with the lowest energy, which is 
not necessarily the ground state especially at the critical 
point with high numerical noise. Afterward, we rotate 
the boson states with a phase angle 0 following the same 
procedure as that adopted for the DMRG calculations. 
This approach is justified because at the critical point 
the system obeys U(l) symmetry. By comparing the en¬ 
ergies of the states, a real ground state can be found sub¬ 
sequently. Details of this approach are provided in the 
Supplementary Material. For each state with 0, we cal¬ 
culate the average displacements X and Z of both baths 
and the value of ( accordingly, as shown in Fig. [4] To a 
large extent, ( is also found to vanish by the variational 
approach, lending support to our DMRG results. More 
importantly, in the deep sub-Ohmic regime, ( completely 
vanishes for a majority of the cases, implying the states 
to be almost orthogonal to each other. The location of 
the real ground state has been indicated by the arrows 
in Fig. in agreement with the DMRG results. More¬ 
over, we have also compared the average displacement of 
all the bosonic modes with the DMRG results and the 
results are presented in Supplementary Material. The 
displacement in the case s = 0.6 is much smaller than 
that of s = 0.25, implying the ground state of s = 0.6 
to be closer to the origin of the X — Z plane rather than 
that of s = 0.25. This finding elucidates why ( is more 
likely to be nonzero for s = 0.6 than for s = 0.25. 
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Conclusion- In summary, we have employed the 
DMRG algorithm to study the TBSBM with both deep 
and shallow sub-Ohmic bosonic spectral densities. The 
numerical approach is adapted with the SOPB, and the 
parity operator is optimally constructed. It is found 
that through this adaptation, the numerical precision 
is greatly improved in both sub-Ohmic regimes. Using 
the SOPB method, we investigate the spin population 
and the ground-state energy. The results show that for 
s = 0.25 both the two quantities change significantly, 
supporting a feature of the first-order phase transition, 
while for s = 0.6 the changes become smoother indicat¬ 
ing a higher order phase transition. We have also put 
forth a newly defined order parameter £, based on which 


various features of the critical phase are discussed and 
the phase diagram is explicitly obtained. It is concluded 
that the SOPB-adapted DMRG algorithm is well-suited 
to handle the complexity of the phase transitions in the 
SBM. This robust approach is expected to be extended 
to tackle other issues, such as the real-time dynamics of 
the SBM. 
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PACS numbers: 


I. OPTIMAL PHONON BASIS 


The calculation procedure of the OPB-adapted DMRG 
algorithm, as sketched in Fig. [Q can be briefly summed 
up as follows. Firstly, we apply the usual DMRG pro¬ 
cedure with restricted basis to obtain the state with the 
lowest energy which may be a rough estimate of the true 
ground state energy of the system (IJ. To accelerate the 
computations, it is useful to add a very small bias e to the 
system in this step to lift the degeneracy. The sign of the 
bias is not important in this step, because in the steps 
later on both degenerate states will be obtained and the 
influence of the small bias will be eliminated by the scal¬ 
ing analysis. Secondly, the DMRG iteration is continued, 
and the basis of the Fock space on the left single site is 
replaced by the bare phonon one during each step of the 
iteration. The number of bare basis ( Nr ) is much larger 
than that of the restricted basis (Nr), ensuring the con- 
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FIG. 1: Schematic for one step of the DMRG iteration with 
OPB. The left single site is initially with bare basis and the 
ground state is calculated for the whole system. Based upon 
the calculated state, the single site with bare basis is opti¬ 
mized. Following that the single site moves forward and the 
above step is redone. The large vertical arrow in the chain 
represents a spin. 
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vergency of the local energy. Thirdly, a new state with a 
much lowered energy with respect to the bare basis is ob¬ 
tained. Based upon this state, we calculate the reduced 
density matrix of the left single site and follow the idea of 
DMRG to discard those bases with low probability. The 
new reduced basis for the local Fock space is the OPB, 
carrying as much information as the bare phonon basis. 
After one iteration, all the boson sites are optimized with 
OPB, and then the system energy is minimized using the 
usual DMRG procedure. Finally, a state with an energy 
very close to that of the real ground state is obtained, 
with which we can calculate all desired observables. 


II. EFFECTIVE DISPLACEMENT OF BOSONS 
A. Deep sub-Ohmic regime (s < 0.5) 


In order to examine the reliability of our numerical ap¬ 
proach, we apply the operator Vi (59) repeatedly onto the 
ground state | g) to obtain a new state |0). As discussed, 
the states with opposite displacements are orthogonal to 
each other, such that the new state with 6 = tt is orthog¬ 
onal to the ground state, i.e., ( g\9 = tt) = 0. Therefore, 
if the information of the state with the opposite displace¬ 
ment is not taken into account in the procedure of the 
numerical optimization, the action of the operator Vi (58) 
will give rise to results completely out of line with ex¬ 
pectations. The examination of this operator is a critical 
criterion to determine whether our numerical approach is 
able to attain the double degeneracy, and thus the phase 
transition. 

Based upon the state |0), we calculate the displace¬ 
ments of every boson mode Xi(8)(= l/\/2(b\ + bf)) with 
both the symmetrically optimized phonon basis (SOPB) 
and asymmetrically optimized phonon basis (AOPB), 
and results shown in Fig. [2] are for the case of a x = 0.02 
and 8 odd sites that are closest to the spin. It is found 
that with SOPB, displacements all switch signs when 
8 = 7r, while the same signs are recovered when 9 = 2n. 
On the other hand, the AOPB displacements are found 
to be smaller than those of SOPB. More importantly, 
after a 2n rotation the displacements fail to return to 
their original positions, implying anomalies in the curves 
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FIG. 2: The effective displacement of boson on the odd 
sites that are closest to the spin versus the rotation angle 
with (a) the SOPB and (b) the AOPB. The parameters are 
s = 0.25, a z = a x = 0.02, and Nb = 16. Different colours 
correspond to different sites, and the arrows indicate the ori¬ 
entation away from the spin. 



0/7C 

FIG. 3: The effective displacement of boson on the odd 
sites that are closest to the spin versus the rotation angle 
with (a) the SOPB and (b) the AOPB. The parameters are 
s = 0.6, a z = a x = 0.1, and Nb = 16. Different colours 
correspond to different sites, and the arrows indicate the ori¬ 
entation away from the spin. 


when 9 becomes large. We also calculate the energy (not 
shown), which is directly related to the displacement of 
the boson modes, and for the two cases we also observe 
that the energy of SOPB case at 9 = 2tt reproduces the 
original value while that of AOPB case does not. These 
results clearly show great improvements when the parity 
symmetry is explicitly taken into account. In addition, 
the deviation |AW| for the displacement of 9 = 27t from 
that of 9 = 0 is a perfect measure of the numerical preci¬ 
sion. The value of |AXi| is found to be smaller than 10 -5 
in Fig.[2j indicating the results in the case of s = 0.25 and 
a x = 0.02 are reliable. We have also checked the value 
of (g\9 = 7r) which is of the same order as the deviation. 


B. Shallow sub-Ohmic regime (0.5 < s < 1) 


In order to indicate the benefits of the SOPB method 
in the shallow sub-Ohmic regime, we show in Fig. [3] the 
phase angle dependence of the displacements for both 
SOPB and AOPB cases similar to that in the previous 
subsection. At variance to the deep sub-Ohmic case, the 
shallow sub-Ohmic bath produces no anomalies for the 
two bases, as the displacements show perfect rotations in 
line with the phase angle implying great precision. On 
the other hand, the AOPB method gives rise to much 
larger displacements than those from the SOPB method. 
Similar with the analysis in the deep sub-Ohmic case, a 
large displacement here would result in a large polariza¬ 
tion of the spin, and the AOPB spin population shows a 
sudden change at the transition point, while the SOPB 
one shows a continuous crossover. The latter case is more 
accurate according to our previous mean-field analysis. 
Therefore, the applicable extent of the method is deter¬ 
mined, that is, the SOPB-adapted method works well 
for both deep and shallow sub-Ohmic cases while dealing 
with the phase transition. 


III. THE VARIATIONAL APPROACH FOR THE 
TBSBM 


A. The rotational optimization 


In order to make a comparison, we adopt the varia¬ 
tional approach [2j to calculate the quantities. The vari¬ 
ational approach is based on the trial wavefunction 
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|0>ph 
|0)phj (1) 


where |+) (|—)) stands for the spin up (down) state in 
the Z-direction, H.c. denotes the Hermitian conjugate, 
|0) p h is the vacuum state of the phonon bath, M and N 
represent the numbers of the bath modes and coherent 
superposition states, respectively, and A n , B n , f n j and 
g U: i are the variational parameters. Herein, the variables 
n and l represent the rank of the coherent superposition 
state and the label of the bosonic mode, respectively. By 
defining the norm of the wave function as D = (4’\4’) and 
the energy E as E = H/D with H = {ijj\H\'i/j), one can 
derive the self-consistency equations as the variational 
conditions, namely 


dH dD 

3 {&} 0 {&} 


( 2 ) 


where denotes one of the variational parameters. The 
ground state is then searched by using this variational 
iteration technology. 

Due to the numerical errors, however, the state | ip) 
found by the variational approach may be a meta-stable 
state. This situation is quite similar to that in the DMRG 
calculation. In particular, since the coherent states are 
used as the trial wave function, the calculated critical 
point will always deviate from the theoretical point a x = 
a z . Motivated by the SOPB-adapted DMRG algorithm, 
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FIG. 4: |(ct z )| calculated by the variational approach as a 
function of a x ranging from 0.018 to 0.021 for s = 0.25 and 
a z = 0.02 with and without the rotational optimization (RO). 
The blue dash-dotted line indicates the critical point at a x = 

Otz- 

we then incorporate the symmetry into the variational 
approach. We construct a rotational operator T(0) = 
exp(i0S'), where 

S = \°y + i {h* b lz - b U b i,*) ( 3 ) 

is the generator of the U(l) symmetry. It is worth noting 
that, there is a hopping term between the two baths in 
the generator S, making it hard to implement the U(l) 
symmetry in the DMRG algorithm which will be dis¬ 
cussed later on. By applying the operator T(0) onto the 
state \ip), we are able to find the ground state with the 
lowest energy in the X — Z plane, which is much more re¬ 
liable than that found by the variational method without 
the optimization. 

In the case with a z = 0.02 and s = 0.25, a sharp jump 
of the magnetization |(<r z )| has been found at the critical 
point a XtC ~ 0.02, similar to the results obtained by the 
DMRG algorithm with SOPB. To locate the critical point 
explicitly, the behavior of the magnetization in a much 
narrower range [0.019, 0.021] is investigated. As shown 
in Fig. [2 the value of the critical point a XyC = 0.02000(1) 
is located with the rotational optimization, an improve¬ 
ment over a XjC = 0.02008(1) obtained without the opti¬ 
mization, since the theoretical value of a x ^ c = a z = 0.02 
is expected. It indicates that the rotational optimization 
makes the variational approach more precise in determin¬ 
ing the critical point. Also calculated is the ground state 
energy E g: as shown in Fig. [5] The slope of the curve 
E g (a x ) is comparable with that obtained by DMRG, 
though the explicit value of E g is slightly larger than 
that from DMRG, e.g., E g = —0.12917 for the variational 
method and E g = —0.12923 for the DMRG algorithm at 
s = 0.6 and a x = a z = 0.1, implying the high precision 
of the DMRG algorithm. Moreover, two distinct values 
of the slopes, 0.89 and 1.55, are obtained for the localized 
and delocalized phases, respectively, indicating a first or¬ 
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X 

FIG. 5: The ground state energy E g as a function of a x at 
s = 0.25 and a z = 0.02 calculated by the variational approach 
with rotational optimization. The dashed lines represent the 
linear fits in both the localized (red) and delocalized phases 
(blue). 


der phase transition. 


B. Order parameter calculated by the variational 
approach with rotational optimization 

The order parameter £ is also calculated by the varia¬ 
tional approach with the rotational optimization, in com¬ 
parison with that from DMRG algorithm with SOPB. 
As shown in Fig. [6] ( is plotted for the two cases of 
s = 0.25, a z = 0.02 and s = 0.6, a z = 0.1. In the 
deep sub-Ohmic regime with s = 0.25 < 0.5, £ is found 
to be vanishing at the critical point ct x /a z = 1.0, i.e, 
£(1.0) = 0. On the other hand, a slight decrease from 
C =1 is found at the critical point in the shallow sub- 
Ohmic regime with s = 0.6. The result is in good agree¬ 
ment with that obtained by DMRG, indicating the po¬ 
tency of the SOPB-adapted DMRG algorithm. 

To further understand the order parameter ( at the 
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FIG. 6: The order parameter £ calculated by the variational 
approach with the rotational optimization as a function of 
the ratio a x /a z for the two cases of s = 0.25 ,a z = 0.02 and 
s = 0.6, a z — 0.1, respectively. 
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FIG. 7: The order parameter £ calculated by the variational 
approach with the rotational optimization is displayed as a 
function of 0/-7T for the two cases of s = 0.25, a x = a z = 0.02 
and s = 0.6, Ox — a z = 0.1 in the deep and shallow sub- 
Ohmic regimes, respectively. The arrows indicate the loca¬ 
tions of the calculated ground states. 


critical point, the behavior of ((0) is investigated on 
the states T(@) \ip) with the rotational optimization, as 
shown in Fig. [ 7 ] For s = 0.25 in the localized phase, 
sharp peaks of £(0) are found at 0/tt « n/2 (n = 0,1, 2) 
with a small peak width A0i defined as the size of the 
parity-symmetry regime with £ > 0. It suggests that the 
ground state is localized in a corner of the X — Z plane, 
where X and Z corresponding to the diagonal and off- 
diagonal coupling baths, respectively. For s = 0.6 > 0.5 
in the critical phases, the width of the peak A 02 in¬ 
creases visibly. Interestingly, the ratio of these two peak 
widths A0 2 /A0 i = 6.0 is consistent with the ratio of 
ground state central angles being 6.3 in the X — Z plane, 
indicating that the width of the peak is determined by 
the central angle. Since the ground state marked by the 
arrow is located nearby the peak, the value of £ in the 
ground state increases with the width of the peak. It 
elucidates why £ is more likely to be nonzero for s = 0.6 
than that for s = 0.25 in Fig. [6] 


IV. U(l) SYMMETRY OF TBSBM WITH 
SINGLE MODE 

At a x = a z , the entire system obeys the U(l) sym¬ 
metry with the generator expressed in Eq. ©■ Very 
recently, Bruognolo et al. tried to incorporate the U(l) 
symmetry explicitly into the DMRG algorithm [3j. It is 
not easy done since there is a hopping term between the 
two bosonic baths in the generator. To implement the 
symmetry, they make a transformation in the Hamilto¬ 
nian of TBSBM to mix the boson modes in the two baths. 
It is subsequently found that the approach does not work 
so well for the case of s < 0.5 due to the large truncation 
error. In comparison, our approach with SOPB can be 
applied to all cases, especially for the cases with s < 0.5 
and a z ^ a x . 


It was stated that in the numerical calculations the 
U(l) symmetry will always be broken spontaneously Q. 
Indeed, our DMRG results do not show simultaneous 
magnetization disappearance in both Z and X orienta¬ 
tions, i.e., (a z ) = (cr x ) = 0. Whereas, the results of £ 
show that the symmetry is still held. This shows that the 
involvement of the bosonic component in the order pa¬ 
rameter is essential. To confirm this statement, we adopt 
the variational approach to investigate the system with 
two single-mode bosonic baths (i.e., M = 1 in Eq. ©), 
which could be dealt with in a numerically exact man¬ 
ner. We will show that in the exact case, the symmetry 
is spontaneously broken in the strong coupling case with 
A » oj. 

The ground state wave function of the TBSBM is also 
investigated to reveal the nature of the localized and de¬ 
localized phases. The phonon population P(x, z) is in¬ 
troduced to describe the ground state of the two-bath 
model. Assuming that the wave function of the ground 
state can be written as 

l^s) = l+)IV’+)ph + |~)|^-)ph, (4) 

where |•0_ ( _) p h and |-0_) p h are the phonon parts of the 
wave function corresponding to the spin up and down 
states, respectively, which can be expanded in a series of 
Fock states or coherent states. Thus the phonon popula¬ 
tion P(x , z) is defined as 

P(x,z) = (ar, z)\'Sf g (x, z))sp in 

= \i>+(x,z)\ 2 + \tp-(x,z)\ 2 (5) 

where (• • • ) sp i n represents the trace of the spin freedom 
of the wave function, x and z correspond to the off- 
diagonal and diagonal coupling baths, respectively, and 
i/j±(x,z) = (f|^±) p h is the phonon part of the wave 
function in the two-dimensional coordinate representa¬ 
tion r = (x, z ). 

According to the definitions of the phonon population 
P(x,z) in Eq. ([5]) and the variational ansatz in Eq. ([Tjl. 
P(x, z) can be calculated as 

P(x z) = y t [AnFnjx, z)] 2 + [B n G n (x,z)] 2 ^ 

n =1 

where D = (’F g |'I' g ) is the norm of the wave function, 
A n and B n denote the weight coefficients of the n-th co¬ 
herent superposition state coupled to the spin up and 
down states, respectively, and F n (x,z) = (r|'i/’+)ph = 
fn,x{x)fn,z{z ) and G n (x,z) = (r|V>_)ph = gn,x{x)g n ,z{z ) 
represent the phonon part of the wave function IV'ilph 
in the coordinate representation f = (x, z). The function 
fn,x{x) denoting the phonon state in the off-diagonal cou¬ 
pling bath then can be obtained, 

fn,x ( X ) = IIW"’*) ( 7 ) 

l 

— e ~ixipi/2 e ipix e ~uj(x-xi) 2 /2, 

l 77 
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In the same way, the functions f n ,z[z),g n , x {x) and g n ,z{z) 
can also be calculated from the displacement coefficients 
f n j and g Ut i. Without loss of generality, we set the fre¬ 
quency uji to unity such that X and Z are now dimen¬ 
sionless. In the single-mode case with the number of 
effective bath modes M = 1, the function F n (x,z) = 
fn,xi x )fn,z{z) = (z\fn,i)(x\fn, 2 ) is simplified, where the 
subscript 1 and 2 correspond to the diagonal and off- 
diagonal coupling baths, respectively. 


FIG. 8: The wave function of the ground state for the single¬ 
mode case (i.e., M=1 in Eq. ©) is plotted in (a) for the 
strong coupling regime with X/uj = 10 and in (b) for the 
weak coupling regime with X/w = 0.1. The X-coordinate 
and Z-coordinate correspond to the off-diagonal and diagonal 
coupling baths, respectively, and the colour represents the 
distribution of the phonon population P(x,z). 


where Xi and pi are defined as 

Pi = ~ i \f^ (fn,l ~ fn,l) > (8) 

Xi = -±= (f + /* ) . (9) 

V 2 w; 


As shown in Fig. [8] (a), the distribution of the phonon 
probability density P(x,z) is located at a corner of the 
X-Z plane for the strong coupling regime with X/uj = 10, 
corresponding to the localized phase. Quite different 
from that in (a), the probability density distribution of 
the phonon P(x, z) in the critical phase concentrates at 
the origin, as shown in Fig. [8] (b) for the weak coupling 
regime with X/ui = 0.1. It shows that the distance be¬ 
tween the center of P(x, z ) and the origin is an important 
parameter to watch as one tells the localize phase from 
the critical one. 
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